Stokes drift and its discontents

The Stokes velocity uS, defined approximately by Stokes (1847, Trans. Camb. Philos. Soc., 8, 441–455.), and exactly via the Generalized Lagrangian Mean, is divergent even in an incompressible fluid. We show that the Stokes velocity can be naturally decomposed into a solenoidal component, usolS, and a remainder that is small for waves with slowly varying amplitudes. We further show that usolS arises as the sole Stokes velocity when the Lagrangian mean flow is suitably redefined to ensure its exact incompressibility. The construction is an application of Soward & Roberts’s glm theory (2010, J. Fluid Mech., 661, 45–72. (doi:10.1017/S0022112010002867)) which we specialize to surface gravity waves and implement effectively using a Lie series expansion. We further show that the corresponding Lagrangian-mean momentum equation is formally identical to the Craik–Leibovich (CL) equation with usolS replacing uS, and we discuss the form of the Stokes pumping associated with both uS and usolS. This article is part of the theme issue ‘Mathematical problems in physical fluid dynamics (part 1)’.


Introduction
Surface gravity waves induce a rectified motion of fluid particles and thus a wave-averaged difference between the mean Eulerian velocity, u E , and the mean Lagrangian velocity u L [1,2] u L = u E + u S .  vortex force and the Stokes-Coriolis force, in the wave-averaged momentum and vorticity equations [3][4][5][6][7].
The Stokes velocity can be defined exactly at finite wave amplitude using Generalized Lagrangian Mean (GLM) theory [8,9]. This exact GLM u S is rotational and compressible, even if the underlying fluid motion is irrotational and incompressible [10]. Expansion in powers of a wave-amplitude parameter produces the standard approximation [2,11] to the Stokes velocity u S = (ξ 1 · ∇) u 1 . (1. 2) The overbar in (1.2) denotes a running time mean, or phase average. In (1.2), u 1 is the linear (first order in ) velocity of the wave and the associated displacement ξ 1 is defined by ∂ t ξ 1 = u 1 and ξ 1 = 0. (The subscript 1 indicates the first-order fields throughout.) The small-amplitude approximation to u S in (1.2) is also rotational and compressible: assuming only that ∇ · ξ 1 = 0, McIntyre [10] shows from (1.2) that The time derivative of an averaged quadratic quantity in (1.3) entails the same slow-modulation assumption that underlies the concept of group velocity and so introduces a second small parameter, μ. The Eulerian mean velocity, u E in (1.1), is incompressible and thus the divergent u S in (1.3) implies a divergent Lagrangian mean velocity.
In figure 1, we illustrate the role of the two small parameters and μ, by considering a weakly nonlinear, slowly modulated two-dimensional packet of deep-water surface gravity waves. The Stokes expansion [1,10] is justified by the weakly nonlinear assumption that the wave slope is small: = ak 1, where k is the wavenumber and a is the amplitude of the surface displacement. In this example the slow-modulation parameter is μ = (k ) −1 1 where is length scale of the packet envelope. Figure 2 shows the motion of a fluid particle in the velocity field of this wave.
Despite (1.3), and the exact results provided by GLM, some authors are reluctant to accept the reality of non-zero ∇ · u S . Moreover, discontent with ∇ · u S = 0 is sometimes confounded with unease over the vertical component of the Stokes velocity, w S =ẑ · u S . For example, rather than taking the vertical component of (1.2), McWilliams et al. [12] defines a 'vertical Stokes pseudovelocity' which, together with the horizontal components of u S , makes an incompressible threedimensional 'Stokes pseudo-velocity'.
Mellor [13], while emphasizing that u S is divergent, is unwilling to accept a non-zero vertical component w S : 'a mean vertical drift is not acceptable'. In the view of concerns with the mean   Panels (a,b) show the x-and z-displacements as functions of time; panel (c) shows the trajectory. In this computation, we assume that the depth d is much greater than the packet length scale so that the second-order Eulerian mean flow is negligible in the wave-active zone. (Online version in colour.) vertical drift w S it is reassuring that particle tracking velocimetry can be used to observe vertical Lagrangian displacements, w S = 0, beneath groups of deep-water waves [14,15]: the mean vertical drift is upward as a wave packet arrives and downwards as the packet departs. (For the wave packet in figure 1, the maximum vertical displacement resulting from w S is about 4.8 cm-this is not visible in figure 2.) After the passage of the packet a fluid particle returns to its initial depth. It is a net vertical displacement that is unacceptable: transient vertical motion, on time scales longer than a 10-s wave period, and shorter than the 100-s packet transit time, is not a concern.
Mellor critiques the Craik & Leibovich ([3,4]) vortex-force formulation of wave-mean interaction by arguing that CL and subsequent authors incorrectly assume that the divergence of u S is zero. A different interpretation is that CL and many followers assume that the wave field has no temporal modulation so that the right of (1.3) is conveniently zero. For example, McWilliams & Restrepo [7] claim to prove ∇ · u S = 0. But examination of this argument shows that [7] assumes that there is no temporal modulation of the wave field. This raises the issue of whether the CL formulation is incomplete or misleading in situations with temporal modulation of the wave field, e.g. in ocean observations [16] and in modelling the growth of swell [17].
In this paper, we revisit the concept of Stokes velocity. For surface gravity waves, we exhibit a natural Helmholtz decomposition of u S in (1.2) and argue that the solenoidal component, u S sol , can advantageously replace u S in most situations. We emphasize that the familiar form (1.2) of the Stokes velocity is not unique but depends on a specific definition of the Lagrangian-mean flow, royalsocietypublishing.org/journal/rsta Phil. Trans. R. Soc. A 380: namely the GLM definition of Andrews & McIntyre [8]. An alternative definition, proposed by Soward & Roberts [18] and closely related to classical averaging and its Lie series implementation (e.g. [19,20]), leads to a solenoidal Lagrangian-mean velocity with u S sol as the corresponding Stokes velocity. This alternative definition, known as 'glm' but better characterized as 'solenoidal Lagrangian mean', has the added benefit of coordinate independence in any geometry, unlike standard GLM (see [21] for other coordinate-independent definitions of the Lagrangian mean). We show that, for surface gravity waves, the associated Lagrangian-mean momentum equation governing the dynamics of the Eulerian mean flow is the CL equation with u S sol replacing u S . The difference between GLM and glm is vividly illustrated with an example proposed by O. Bühler (personal communication, 2021). Consider a bucket of initially motionless water. If the water is agitated, for example, by pressure forcing at the surface, then the potential energy of the water is increased, or equivalently the centre of mass of the water is elevated above its initial height. The fluid at the bottom of Bühler's bucket, however, cannot move in the vertical and so any definition of 'Lagrangian mean' that tracks the position of the centre of mass of the fluidsuch as GLM-will be divergent in this situation, even though the velocity of the water is entirely incompressible. Conversely, any definition of 'Lagrangian mean' resulting in a strictly solenoidal Lagrangian mean velocity-such as glm-cannot track the centre of mass of the fluid.
The plan of the paper is as follows. In §2, we sketch the derivation of the standard form (1.2) of the Stokes velocity, give its Helmholtz decomposition, and show how a simple modification of this derivation, implementing an alternative Lagrangian-mean flow definition, naturally brings about the velocity u S sol . In §3, we examine the respective role of u S and u S sol in Stokes pumping, which is the mechanism whereby the horizontal divergence of the Stokes transport drives an Eulerian mean flow. In §4, we show how the glm approach enables the systematic construction of solenoidal Lagrangian-mean and Stokes velocities up to arbitrary algebraic accuracy in . We explain how Lie series provide both an interpretation and an efficient implementation of this construction, and we derive the glm version of the CL equations. Section 5 gives the conclusion.

The Stokes velocity u S and its solenoidal part u S sol (a) Derivation of the Stokes velocity
We start by recalling the traditional derivation of the Stokes velocity in (1.2). The position x(t) of a fluid particle is determined by solving where α is a wave phase, regarded as an ensemble parameter. The fluid velocity u(x, t, α, ) is incompressible, ∇ · u = 0, and has the form where is the wave amplitude parameter. The leading-order term u 1 (x, t, α) is a fast wavy flow, so u 1 = 0, where the mean, denoted by the overbar, is an average over the phase α.
To average the fast wave oscillations in (2.1), we consider the ansatz is the slow motion of a Lagrangian mean position. Think of x L as a 'guiding centre' such that rapid wavy oscillations are confined to the displacements ξ n ; these displacements from x L do not grow with time, i.e. all members of the ensemble remain close to the guiding centre x L . The motion of x L is written as where u L (x, t, ) is the Lagrangian mean velocity, yet to be defined and determined. The ansatz (2.3b) is ambiguous because requiring only that the ξ n 's do not grow with time does not uniquely determine u L and ξ n . We return to this point below. Substituting (2.3b) and (2.5) into (2.1) and (2.2) and matching powers of at the first two orders results in Choosing to follow Stokes [1] and Andrews & McIntyre [8], one disambiguates (2.3b) by requiring that ξ 2 (x L , t, α) = 0. In this case, averaging (2.6b) produces the familiar result In (2.7) we have now omitted the subscript 2 on u L . With (2.7) we recover (1.2) and the small wave-amplitude version of (1.1).

(b) The solenoidal Stokes velocity u S sol
Let us examine the divergence of u S in more detail. The 'un-averaged Stokes velocity' can be written exactly as In passing from (2.8a) to (2.8b), we have used wave incompressibility to simplify the standard vector identity for the curl of the cross product u 1 ×ξ 1 . The average of (2.8b), identifies a solenoidal part of the Stokes velocity as The solenoidal vector u S sol is the incompressible part of the Stokes velocity for all types of weakly nonlinear waves in an incompressible fluid. We propose that u S sol can advantageously replace the traditional form of the Stokes velocity (1.2) in many circumstances.
The solenoidal Stokes velocity arises naturally if a small change is made in the derivation in §2(a): suppose we decide from the outset to change the definition of 'Lagrangian mean' so that (rather than ξ 2 = 0). Averaging (2.6b) then results in The incompressible Lagrangian mean velocity in (2.12) is an alternative to the traditional compressible u L in (2.7). The O( 2 ) shift ξ 2 in the position of guiding centre x L (t, ) produces Lagrangian mean and Stokes velocities that are divergence free at order 2 . In §4, we discuss a systematic framework-Soward & Robert's [18] glm alternative to GLM-that generalizes this property to arbitrarily high order in . Note that the difference between u S and u S sol is a time derivative so has no impact on particle dispersion for waves that are represented by stationary random processes as discussed by Holmes-Cerfon & Bühler [22].

Stokes transport and Stokes pumping
Discussions of the deep return flow associated with a surface-gravity wave packet argue [2] that: 'the return flow · · · can be explained as the irrotational response to balance the Stokes transport · · · that acts to "pump" fluid from the trailing edge of the packet to the leading edge'. Similar sentiments are expressed in [23,24]. With this physical picture in mind, it is instructive to compare the Stokes pumping associated with u S with that of u S sol . The interpretation of the return flow in the quotation above is most consistent with u S sol . (a) Stokes pumping and u S sol = ∇× 1 2 u 1 × ξ 1 We start with the easy solenoidal case, with Stokes transport where the vertical integration above is from the bottom of the wave-active zone (denoted −∞) to the mean sea surface at z = 0. (In this section, we confine attention to deep-water waves so that the lower limit, −∞, is well above the distant bottom.) Vertical integration of ∇ · u S sol = 0 over the wave-active region produces the unsurprising result where we use 0 to denote evaluation at z = 0, e.g. w S sol 0 = w S sol (x, y, 0, t). The result (3.2b) is consistent with the idea that horizontal convergence within the wave-active zone pumps fluid downwards, out of the wave-active zone, with the vertical velocity w S sol (x, y, 0, t). Vertical integration of u S sol = 1 2 (ξ 1 · ∇)u 1 − 1 2 (u 1 · ∇)ξ 1 over the wave-active zone results in is a streamfunction for horizontally circulating Stokes transport. This horizontal Stokes circulation is previously unremarked, perhaps because the terms involving χ in (3.3a) and (3.3b) are order μ smaller than the other terms. Using the coordinate expressions in (3.3a) and (3.3b), the solenoidal pumping can be expressed entirely in terms of surface quantities: (b) Stokes pumping and u S = (ξ 1 · ∇)u 1 We turn now to the traditional definition of the Stokes velocity. Define the Stokes transport via Because ∇ · u S = 0 there is no analogy to (3.2b). Instead, integrating (2.15) over the wave-active zone With vertical integration of the horizontal components of u S over the wave-active zonê (3.10) In the final term, the indices α and β are 1 and 2. Eliminating ∇ · T S between (3.8) and (3.10) The results in (3.8), (3.10) and (3.11) are all more complicated than their solenoidal cousins in (3.2b), (3.3a) and (3.3b).

(c) A comment on approximations to the Stokes transport
A widely used expression for the Stokes transport is Expression (3.12) is exact for a uniform (μ = 0) progressive wave and is a leading-order approximation for a slowly modulated (μ 1) wave packet. To see the connection with the more general and exact results above, align the x-axis with the horizontal wavevector so that η 1 = v 1 = 0; thenŷ · T S =ŷ · T S sol = χ = 0. Thus, with an error by order 2 μ, and in agreement with (3.12), In other words, (3.12) is a leading-order approximation and at this order T S and T S sol are identical.

glm
We return to the definition of the Stokes velocity and show how the glm theory of Soward & Roberts [18] rationalizes and generalizes the heuristic construction of the solenoidal velocity u S sol . In §2(a), we emphasized the ambiguity in the decomposition (2.3a) of trajectories into a mean part x L and a perturbation ξ . GLM resolves this ambiguity by imposing that While (4.1) is widely accepted, it is neither inevitable nor particularly natural. It implies that the Lagrangian-mean trajectory of a particle is defined by the equality x L = x between the coordinates of the Lagrangian-mean position and the average of the coordinates of the particle. This indicates that the Lagrangian-mean trajectory and, as a result, the Lagrangian-mean velocity u L depend on a choice of coordinates. This undesirable feature of GLM is remedied by glm, while also ensuring that u L is non-divergent (see [21] for other alternatives to GLM).

(a) Formulation
To introduce glm, it is convenient to rewrite equation (2.1) governing the fluid trajectories in terms of the flow map ϕ(x, t, α, ) giving the position at time t of the fluid particle initially at x.  t, α, ) = u(ϕ(x, t, α, ), t, α, ). (4.2) The decomposition of trajectories into mean and perturbation is best written as the composition ϕ = Ξ • ϕ L , (4.3) or, more explicitly, ϕ(x, t, α, ) = Ξ (ϕ L (x, t, ), t, α, ). Here ϕ L is the (α-independent) mean map, sending the initial position of particles to their mean position, and Ξ is the perturbation map, sending the mean position to the exact, perturbed position. The perturbation map can only be represented in the familiar form x → x + ξ (x, t, α, ), as in (2.3a), in Euclidean space, where positions x can be identified with vectors and added, or, in more general geometries, once a specific coordinate system has been chosen and x is interpreted as a triple of coordinates. The smallness of the perturbation, usually stated as |ξ | 1, translates into the requirement that Ξ is close to the identity map. We emphasize that the mean map ϕ L is not obtained from ϕ by applying an averaging operator: averaging is a linear operation that applies to linear objects, such as vector fields, but not to nonlinear maps such as ϕ. Instead, ϕ L is defined by imposing a condition on the perturbation map Ξ . The form of the Lagrangian-mean velocity u L , defined bẏ depends on this condition. The defining condition of glm is expressed as follows. The small parameter is regarded as a fictitious time, and the perturbation map Ξ is constructed as the flow at 'time' of a vector field, q say; that is, (4.5) and t is treated as a fixed parameter. The glm condition is then Although superficially similar to the GLM condition (4.1), (4.6) is fundamentally different in that it is an intrinsic statement, applicable to any manifold and independent of any coordinate choice. Moreover, the glm formulation defines an exactly divergence-free Lagrangian mean flow, by requiring that ∇ · q = 0 (4.7) to ensure that Ξ and ϕ L preserve volume and hence ∇ · u L = 0. Equation (4.6) generalizes the condition on ξ 2 in (2.12) which leads to u S sol as the Stokes velocity. We show this by solving (4.5) by Taylor expansion. In coordinates, we have that (4.8) The power series expansions (2.3b) for ξ and q = q 1 + q 2 + 2 q 3 + · · · (4.9) for q then give ξ 1 = q 1 and ξ 2 = q 2 + 1 2 q 1 · ∇q 1 (4.10) so that (4.6) implies (2.12). In the next section, we develop a systematic computation of u L and hence u S sol order by order in using Lie series. This removes the need to introduce ξ by focusing the perturbation expansion on q.

(b) Lie series expansion
The glm formalism can be regarded as an instance of classical perturbation theory, which approximates solutions to the ordinary differential equation the variable transformation and ϕ L represents the new variable. Lie series [19,20] provide a powerful tool for the systematic implementation of classical perturbation theory which we now apply to glm.
Introducing the decomposition (4.3) into (4.2) and using (4.4) gives w + Ξ * u L = u, (4.11) where is the perturbation velocity and Ξ * is the push-forward by Ξ , with Ξ * u L = (u L · ∇)Ξ in Cartesian coordinates. Pulling back (4.11) gives We seek an -dependent Ξ to eliminate fast time dependence from u L order-by-order in , and formulate the problem in terms of the vector field q that generates Ξ according to (4.5). We impose that ∇ · q = 0 to ensure that ∇ · u L = 0, and the glm condition (4.6). Expanding q as in (4.9) we relate the various terms by differentiating (4.13) repeatedly with respect to and evaluating the results at = 0. Two key identities turn this into a mechanical exercise. The first, essentially the definition of the Lie derivative [25], is (4.14) where L q u = q · ∇u − u · ∇q is the Lie derivative of u along q. The second, is established in appendix A. Iterating (4.14) and (4.15) we find Introducing (4.16) and (4.17) into (4.13) gives at the first two orders in , ∂ t q 1 = u 1 , i.e. q 1 = ξ 1 (4.18) and on using (4.6). This provides the geometric formula for the solenoidal Stokes drift, equivalent to (2.11a) since L q 1 u 1 = q 1 · ∇u 1 − u 1 · ∇q 1 = ∇ × (u 1 × ξ 1 ). The expansion can be pursued to higher orders by choosing divergence-free q n that push the fast dependence on the right-hand side of (4.13) to order n+1 . This yields u L and hence u S sol to arbitrary order in . on the case of small-amplitude surface gravity waves. This derivation is conveniently carried out using a representation of the rotating Euler equation (4.21) with D t = ∂ t + u · ∇, in terms of the absolute momentum 1-form [21,26,27] ν a = u · dx + 1 2 (f × x) · dx. (4.22) It can be checked that (4.21) is equivalent to where π = p − 1 2 |u| 2 − 1 2 (f × x) · u, using basic properties of the Lie derivative, namely Leibniz rule, commutation with the differential d, and that L u = u · ∇ when applied to scalars [25]; see appendix A for details. Equation (4.23) can be thought of as a local version of Kelvin's circulation theorem. This is readily obtained by integration along a closed curve C(t) moving with u to find An advantage of (4.23) is that it leads directly to a Lagrangian-mean momentum equation of a similar form [21], (4.25) and to the corresponding Lagrangian-mean Kelvin's circulation theorem C L (t) ν L a = const, (4.26) where the closed curve C L (t) moves with the Lagrangian-mean velocity u L . Here the Lagrangianmean momentum and effective pressure are given by ν L a = Ξ * ν a and π L = Ξ * π. (4.27) The pull-back Ξ * by the perturbation map acts on scalars as a composition, e.g. (Ξ * π )(x, t, α, ) = π (Ξ (x, t, α, ), t, α, ), and commutes with the differential so that We show in appendix A that (4.25), together with the coordinate representation Ξ (x, t, α, ) = x + ξ (x, t, α, ) and GLM condition (4.1) recovers Andrews & McIntyre's Lagrangian-mean momentum equation [8, theorem I]. For glm, we can use the Lie-series expression (4.16) to obtain an expansion of ν L a as using that u = u 1 + 2 u 2 + · · · and q 1 = q 2 = 0. Introducing (4.29) into (4.25) and expanding the Lie derivative gives an evolution equation for the Eulerian mean velocity u E = u 2 supplemented by the incompressibility condition ∇ · u E = 0. in (4.29) against the other two in the average (recall that ∂ t q 1 = u 1 ). Moreover, using the irrotationality condition (2.13), we compute L q 1 (u 1 · dx) = (q 1 · ∇)u 1 · dx + u 1 · dq 1 = (q 1j u 1i,j + u 1j q 1j,i ) dx i Therefore, to leading order, the glm mean momentum equation reduces to with the glm Lagrangian-mean velocity in (2.12) (see [26,27] for an analogous formulation of the GLM CL equation). Using Cartan's formula in the form (4.32) the Lie derivatives in (4.31) can be written as where we do not detail the exact differentials. This makes it possible to rewrite (4.31) as for a suitable definition of the effective pressure . Equation (4.35) can be recognized as the CL equation [3,4,28] with the solenoidal Stokes velocity u S sol replacing u S . This is not surprising since the O( 2 μ) difference between u S sol and u S is of the same order as terms neglected in the derivation of the CL equation. (CL geared μ to by taking μ = 2 . In our derivation, this gearing is not necessary.) In fact the assumption μ 1 is only used to neglect the term (4.30) from the Lagrangian-mean momentum equation to obtain (4.31) and hence (4.35). Restoring this term leads to the generalization of the CL equation valid for μ = O(1). Since L u L and d commute, L u L d(∂ t |q 1 | 2 ) = d(· · · ) and the additional terms involving |q 1 | 2 can be absorbed in the differential of the pressure-like term on the right-hand side. We conclude that the CL equation (4.35) with u S sol as Stokes velocity holds for μ = O(1) provided that the effective pressure is suitably redefined.

Conclusion
This paper examines a problematic aspect of the Stokes velocity u S , namely its divergence ∇ · u S and non-zero vertical component w S . Some confusion has arisen because ∇ · u S and w S are small when the wave field has an amplitude that varies on scales longer and slower than the wavelength and period. This scale-separation approximation corresponds to the existence of a small parameter μ 1, as is the case for slowly varying wavepackets. The distinction between approximate and exact results in the literature is not always made plain. Beyond this, there are no irretrievable difficulties with the familiar form (1.2) of u S : when ∇ · u S and w S cannot be neglected, they are readily computed in terms of the first-order fields.
The main point of the paper, however, is that the Stokes velocity, understood as the difference u S = u L − u E between Lagrangian-and Eulerian-mean velocities, is not uniquely defined and that an alternative version, u S sol in (2.11), that is exactly solenoidal can serve as a convenient substitute for the familiar (1.2). The non-uniqueness arises because there is no single definition of the Lagrangian-mean velocity u L , which is only constrained to serve as a good 'guiding-centre' Expanding, we find (D t u) · dx + u · du + 1 2 (f × u) · dx + 1 2 (f × x) · du = −∇p · dx + u · du + 1 2 (f × dx) · u + 1 2 (f × x) · du, ( A 5 ) which recovers (4.21) since (f × dx) · u = −(f × u) · dx.